import pathlib
import sys

# make coastpy library importable by appending root to path
cwd = pathlib.Path().resolve()
proj_dir = cwd.parent.parent
sys.path.append(str(proj_dir / "src"))
# definitly not a best practice, but a workaround to avoid many warnings caused by shapely 2.0 release
import io
import socket
import warnings

import cartopy.crs as ccrs
import geopandas as gpd
import geoviews as gv
import geoviews.tile_sources as gts
import geoviews.tile_sources as gvts
import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import panel as pn
import requests
import rioxarray
import shapely
import xarray as xr
from tqdm import tqdm

from coastpy.geometries import geo_bbox

warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)
/var/folders/p8/qxnqn2ns5kbczrp91kx4_1nm0000gn/T/ipykernel_17541/1733808755.py:7: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:

import os
os.environ['USE_PYGEOS'] = '0'
import geopandas

In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
  import geopandas as gpd
DATA_DIR = proj_dir / "data"
coastsat_metadata_pacific_fp = DATA_DIR / "03_CoastSat_metadata_layer_Pacific.geojson"
metadata = gpd.read_file(coastsat_metadata_pacific_fp)

for c in ["Tier1Images", "Tier2Images", "S2Images", "pk"]:
    metadata[c] = metadata[c].astype("i8")

for c in ["TideRange", "WaveHeight", "RelativeTideRange"]:
    metadata[c] = metadata[c].astype("f8")

metadata.head()
Tier1Images Tier2Images S2Images TideRange WaveHeight RelativeTideRange pk geometry
0 410 113 207 1.65 1.22 1.35 0 POINT (180.00000 -16.15293)
1 256 33 214 1.79 0.60 2.97 1 POINT (179.16992 -16.41845)
2 255 56 415 1.86 0.81 2.29 2 POINT (178.51402 -16.78987)
3 301 78 214 1.47 0.60 2.43 3 POINT (179.18860 -16.72214)
4 495 176 214 1.51 1.01 1.49 4 POINT (179.89177 -16.65793)

Part 1: compute relative tidal range with respect to waves#

todo: check with Kilian which waves heights he used mean_tidal_range / mean_wave_height

metadata["RelativeTideRange"] = metadata["TideRange"] / metadata["WaveHeight"]
# metadata["RelativeTideRange"] = ...  # your code to compute tidal range over the wave height
pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Tide versus wave dominated coasts",
        style={"color": "black"},
        width=500,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

w_threshold = pn.widgets.FloatSlider(
    name="Threshold", start=0.1, end=10.0, step=0.1, value=1
)


@pn.depends(w_threshold.param.value)
def plot(threshold):
    wave_dominated = metadata[metadata["RelativeTideRange"] < threshold].copy()
    tide_dominated = metadata[metadata["RelativeTideRange"] > threshold].copy()

    return (
        gts.EsriImagery()
        * wave_dominated.hvplot(geo=True, color="blue", label="Wave dominated")
        * tide_dominated.hvplot(geo=True, color="green", label="Tide dominated")
    ).opts(legend_position="bottom_right")


app = pn.Column(
    title_bar,
    pn.Row(w_threshold),
    pn.Row(plot),
)
app
from ipyleaflet import Map, basemaps

m = Map(basemap=basemaps.Esri.WorldImagery, scroll_wheel_zoom=True)
m.center = 37.768137, -122.511066
m.zoom = 16
m.layout.height = "800px"
m
bbox = geo_bbox(m.west, m.south, m.east, m.north)
bbox.explore()
Make this Notebook Trusted to load map: File -> Trust Notebook
transects = gpd.read_file(DATA_DIR / "03_CoastSat_transect_layer.geojson")
transects_roi = gpd.sjoin(transects, bbox)
transects_roi.explore()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[10], line 1
----> 1 transects_roi.explore()

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/geodataframe.py:1981, in GeoDataFrame.explore(self, *args, **kwargs)
   1978 @doc(_explore)
   1979 def explore(self, *args, **kwargs):
   1980     """Interactive map based on folium/leaflet.js"""
-> 1981     return _explore(self, *args, **kwargs)

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/geopandas/explore.py:366, in _explore(df, column, cmap, color, m, tiles, attr, tooltip, popup, highlight, categorical, legend, scheme, k, vmin, vmax, width, height, categories, classification_kwds, control_scale, marker_type, marker_kwds, style_kwds, highlight_kwds, missing_kwds, tooltip_kwds, popup_kwds, legend_kwds, map_kwds, **kwargs)
    363         map_kwds["max_zoom"] = tiles.get("max_zoom", 18)
    364         tiles = tiles.build_url(scale_factor="{r}")
--> 366 m = folium.Map(
    367     location=location,
    368     control_scale=control_scale,
    369     tiles=tiles,
    370     attr=attr,
    371     width=width,
    372     height=height,
    373     **map_kwds,
    374 )
    376 # fit bounds to get a proper zoom level
    377 if fit:

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/folium/folium.py:274, in Map.__init__(self, location, width, height, left, top, position, tiles, attr, min_zoom, max_zoom, zoom_start, min_lat, max_lat, min_lon, max_lon, max_bounds, crs, control_scale, prefer_canvas, no_touch, disable_3d, png_enabled, zoom_control, **kwargs)
    272     zoom_start = 1
    273 else:
--> 274     self.location = validate_location(location)
    276 Figure().add_child(self)
    278 # Map Size Parameters.

File ~/mambaforge/envs/codebook/lib/python3.10/site-packages/folium/utilities.py:78, in validate_location(location)
     71         raise ValueError(
     72             "Location should consist of two numerical values, "
     73             "but {!r} of type {} is not convertible to float.".format(
     74                 coord, type(coord)
     75             )
     76         )
     77     if math.isnan(float(coord)):
---> 78         raise ValueError("Location values cannot contain NaNs.")
     79 return [float(x) for x in coords]

ValueError: Location values cannot contain NaNs.
TIMEOUT_SEC = 5  # default timeount in seconds

# Get the transect id's by saving the transects in our region of interest as a list.
transect_ids = transects_roi.TransectId.to_list()

# Here we download the data per transect, but if it takes too long we default to Ocean beach, California.
try:
    data = []
    for transect_id in tqdm(transect_ids):
        url = url = f"http://coastsat.wrl.unsw.edu.au/time-series/{transect_id}/"
        # Pandas.read_csv() can handle url's but it doesn't accept a timeout argument,
        # so we download the data via requests.
        r = requests.get(url, timeout=TIMEOUT_SEC).content
        s = pd.read_csv(
            io.StringIO(r.decode("utf-8")),
            header=None,
            names=["date", "shoreline_position"],
        )
        s["date"] = pd.to_datetime(s["date"])
        s = s.set_index("date").replace({"None": np.nan}).astype("f8")
        s.name = transect_id
        data.append(s)
    column_names = [s.name for s in data]
    shorelines = pd.concat(data, axis=1)
    shorelines.columns = column_names

except:
    # If it takes too long to get the data from the server, default to Ocean Beach, California.
    # First set the map to that area, so that we have the correct bbox.
    m.center = 37.768137, -122.511066
    m.zoom = 16
    m.layout.height = "800px"
    bbox = geo_bbox(m.west, m.south, m.east, m.north)
    transects_roi = gpd.sjoin(transects, bbox)
    shorelines_fp = DATA_DIR / "03_shorelines_ocean_beach.parquet"
    shorelines = pd.read_parquet(shorelines_fp)
  0%|                                         | 0/8 [00:00<?, ?it/s]
shorelines.head()
usa_CA_0215-0003 usa_CA_0215-0004 usa_CA_0215-0005 usa_CA_0215-0006 usa_CA_0215-0007 usa_CA_0215-0008 usa_CA_0215-0009 usa_CA_0215-0010 usa_CA_0215-0011 usa_CA_0215-0012 usa_CA_0215-0013 usa_CA_0215-0014 usa_CA_0215-0015 usa_CA_0215-0016 usa_CA_0215-0017 usa_CA_0215-0018
date
1984-05-02 18:13:15 71.4 87.6 103.9 109.0 115.4 NaN 120.1 124.8 NaN NaN NaN NaN 148.7 146.1 131.7 NaN
1984-06-03 18:14:10 74.4 70.1 98.5 116.8 117.9 141.0 127.4 134.9 NaN 166.6 NaN NaN NaN 139.3 133.6 NaN
1984-06-10 18:20:29 62.9 71.4 90.4 118.9 111.7 137.4 140.6 126.9 137.5 158.1 150.3 150.1 155.8 150.8 142.7 NaN
1984-06-26 18:20:42 60.7 74.3 91.6 112.9 127.7 133.8 152.4 156.0 157.7 171.2 164.9 174.4 165.9 148.2 144.7 NaN
1984-09-07 18:16:09 77.3 91.1 110.3 128.9 141.6 152.9 166.0 173.3 173.5 173.8 174.8 184.2 195.9 202.3 190.7 185.3

#

check this for inspiration on how to clor the plot> https://foundations.projectpythia.org/core/xarray/enso-xarray.html#

pn.extension()
title_bar = pn.Row(
    pn.pane.Markdown(
        "## Historical shoreline series",
        style={"color": "black"},
        width=500,
        sizing_mode="fixed",
        margin=(10, 5, 10, 15),
    ),
    pn.Spacer(),
)

options = transects_roi.TransectId.to_list()
transect_slider = pn.widgets.Select(
    name="Transect", options=options, value=np.random.choice(options)
)


@pn.depends(transect_slider.param.value)
def plot_transects(transect_id):
    transect = transects_roi.loc[transects_roi["TransectId"] == transect_id].copy()

    return (
        gts.EsriImagery()
        * transects_roi.hvplot(geo=True, color="blue")
        * transect.hvplot(geo=True, color="red")
    )


@pn.depends(transect_slider.param.value)
def plot_series(transect_id):
    shoreline_selected = shorelines.hvplot(x="date", y=[transect_id], color="red")
    shorelines_ = shorelines[
        shorelines.columns.difference([transect_slider.value])
    ].hvplot(x="date", alpha=0.4)

    return shorelines_ * shoreline_selected


app = pn.Column(
    title_bar,
    pn.Row(transect_slider),
    pn.Row(pn.Column(plot_transects), plot_series),
)
app
mei_fp = DATA_DIR / "03_meiv2.data"
mei = (
    pd.read_csv(mei_fp, delim_whitespace=True, index_col=0, header=None)
    .reset_index()
    .rename({0: "year"}, axis="columns")
    .melt(id_vars="year", var_name="month")
)
mei["date"] = pd.to_datetime(mei["year"].astype(str) + "-" + mei["month"].astype(str))
mei = (
    mei.drop(["year", "month"], axis="columns")
    .set_index("date")
    .sort_index()
    .replace({-999: np.nan})
)
mei.hvplot(x="date")
mei_interp = mei.reindex(shorelines.index, method="nearest")
shorelines_ = shorelines.copy()
shorelines_.index = mei_interp.value.to_list()
shorelines_.hvplot(kind="scatter")
shorelines_ = shorelines.diff().copy()
shorelines_.index = mei_interp.value.to_list()
shorelines_.hvplot(kind="scatter", alpha=0.2, color="blue")